Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/misc/volume4/dos-fft

Text File | 1989-02-03 | 11.7 KB | 505 lines

Path: xanth!mcnc!gatech!bloom-beacon!husc6!m2c!necntc!ncoast!allbery
From: [e-mail] (Steve Sampson)
Newsgroups: comp.sources.misc
Subject: v04i009: An FFT program
Message-ID: <[e-mail]>
Date: 1 Aug 88 01:44:07 GMT
Sender: [e-mail]
Reply-To: [e-mail] (Steve Sampson)
Organization: The Unix(R) Connection, Dallas, Texas
Lines: 492
Approved: [e-mail]

Posting-number: Volume 4, Issue 9
Submitted-by: "Steve Sampson" <[e-mail]>
Archive-name: dos-fft

[MS-DOS C of some kind -- you may have to patch "fopen"'s under Unix.  At
least.  No binaries were provided.  ++bsa]

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of shell archive."
# Contents:  readme.doc gen.c fft.c
# Wrapped by sampson@killer on Sun Jul 31 20:40:12 1988
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f readme.doc -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"readme.doc\"
else
echo shar: Extracting \"readme.doc\" \(3814 characters\)
sed "s/^X//" >readme.doc <<'END_OF_readme.doc'
XI originally saw an FFT program in Byte Magazine many years ago.  I wrote
Xa version for BASIC that worked pretty good.  Then I thought I'd translate
Xit into C.  These programs are the result.  I don't do windows though...
X
XNeeds a graphic interface, but the time escapes me.  Please upload any better
Xgraphic versions.
X
XThe original Byte Magazine program was designed for real data only.  In my
Xexperiments I needed to preserve both real and imaginary data.  If you feed
Xthe FFT real data only, then the output will be a mirror image, and you can
Xignore the left side.
X
XFor an example try:
X
X    gen 16 in
X    1000
X    3000
X
XWhich will sample the 1 Khz data every 333 microseconds (1 / 3 Khz).
XNote: The sample frequency should be greater than 2 times the input
Xfrequency (Nyquist and all that...).
X
XThen run fft.exe like so:
X
X    fft 16 in out
X
XAnd you should see a display like so:
X
X0    |=======                      (-1500.0 Hz)
X1    |=====                          (-1312.5 Hz)
X2    |====                          (-1125.0 Hz)
X3    |====                           (-937.0 Hz)
X4    |===                           (-750.0 Hz)
X5    |===                           (-562.5 Hz)
X6    |===                           (-375.0 Hz)
X7    |===                           (-187.5 Hz)
X8    |====        <-------   DC            (000.0 Hz)
X9    |====        <-------   Fundamental        (187.5 Hz)
X10    |======        <-------   Second Harmonic    (375.0 Hz)
X11    |========                    (562.5 Hz)
X12    |==============                    (750.0 Hz)
X13    |========================================================
X14    |============================            (1125.0 Hz)    ^
X15    |===========                    (1312.5 Hz)    |
X                                |
X                [13 - 8 (center)] * 187.5 = 937.0 Hz
X
XThe fundamental display frequency is:
X
X    T  = Time Increment Between Samples
X    N  = Number Of Samples
X    Tp = N * T
X
X    Then F = 1 / Tp
X
X    In the example above, the time increment between samples is
X    1 / 3000 or 333 microseconds.  N = 16, so Tp = 5333 microseconds
X    and 1 / .005333 is 187.5 Hz.
X
X    Therefore each filter is a multiple of 187.5 Hertz.  Filter 8 in this
X    example is center, so that would be zero, 9 would be one, etc.
X
XIn this case the sample interval didn't work so good for the frequency and
Xshows the limitation of the Discrete Fourier Transform in representing a
Xcontinuous signal.  A better sample rate for 1000 Hz would be 4000 Hz,
Xin which case T = 250 ms, Tp = 4 ms, and F = 250 Hz.  1000 / 250 = 4.  The
Xpower should all be in filter 12 (8 + 4) in this case.
X
XLet's run it and see:
X
X    gen 16 in
X    1000
X    4000
X
X    fft 16 in out
X
X0    |
X1    |
X2    |
X3    |
X4    |
X5    |
X6    |
X7    |
X8    |
X9    |
X10    |
X11    |
X12    |========================================================
X13    |
X14    |
X15    |
X
XWell what do you know...
X
XThe output file data can be used by other programs as needed.
X
XBy using negative frequencies in gen.exe you can generate opening targets:
X
X    gen 16 in
X    -1000
X    3000
X    fft 16 in out
X
XProduces:
X
X0    |=======
X1    |===========
X2    |============================
X3    |=======================================================
X4    |==============
X5    |========
X6    |======
X7    |====
X8    |====        <-------- Zero Hertz (DC)
X9    |===
X10    |===
X11    |===
X12    |===
X13    |====
X14    |====
X15    |=====
X
XYou can see in these examples where weighting functions would be nice.
XFor an example of what happens when the imaginary data is not input
X(ie, zeros put in) for a 1000 Hz frequency at 3000 Hz sample rate:
X
X0    |===============
X1    |==================
X2    |===================================
X3    |========================================================
X4    |===========
X5    |====
X6    |==
X7    |=                Trash this part
X---------------------------------------------------------------------
X8    |
X9    |=
X10    |==
X11    |====
X12    |===========
X13    |=======================================================
X14    |===================================
X15    |==================
X
XThe left side is redundant and can be deleted.  This is what the original
XByte Magazine article did (December 1978 Issue).
X
XGood luck, have fun with it,
XSteve Sampson, CompuServe: 75136,626  Unix: [e-mail]
END_OF_readme.doc
if test 3814 -ne `wc -c <readme.doc`; then
    echo shar: \"readme.doc\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f gen.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"gen.c\"
else
echo shar: Extracting \"gen.c\" \(1292 characters\)
sed "s/^X//" >gen.c <<'END_OF_gen.c'
X/*
X *    gen.c
X *
X *    C Version 1.0 by Steve Sampson, Public Domain
X *
X *    This program is used to generate time domain sinewave data
X *    for fft.c.  If you want an opening target - negate the test frequency
X *
X *    Usage: gen samples output
X */
X
X#include <stdio.h>
X#include <alloc.h>
X#include <math.h>
X
X#define    PI2    ((double)2.0 * M_PI)
X
Xmain(argc, argv)
Xint    argc;
Xchar    *argv[];
X{
X    FILE    *fp;
X    double    sample, freq, time, *real, *imag;
X    int    loop, samples;
X
X    if (argc != 3)  {
X        printf("Usage: gen samples output_file\n");
X        printf("Where samples is a power of 2\n");
X        exit(-1);
X    }
X
X    if ((fp = fopen(argv[2], "wb")) == (FILE *)NULL)  {
X        printf("Unable to create write file\n");
X        exit(-1);
X    }
X
X    samples = abs(atoi(argv[1]));
X
X    real = (double *)malloc(samples * sizeof(double));
X    imag = (double *)malloc(samples * sizeof(double));
X
X    printf("Input The Test Frequency (Hz) ? ");
X    scanf("%lf", &freq);
X    printf("Input The Sampling Frequency (Hz) ? ");
X    scanf("%lf", &sample);
X    sample = (double)1.0 / sample;
X
X    time = (double)0.0;
X    for (loop = 0; loop < samples; loop++)  {
X        real[loop] =  sin(PI2 * freq * time);
X        imag[loop] = -cos(PI2 * freq * time);
X        time += sample;
X    }
X
X    fwrite(real, sizeof(double), samples, fp);
X    fwrite(imag, sizeof(double), samples, fp);
X
X    fclose(fp);
X    putchar('\n');
X}
X
X/* EOF */
END_OF_gen.c
if test 1292 -ne `wc -c <gen.c`; then
    echo shar: \"gen.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f fft.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"fft.c\"
else
echo shar: Extracting \"fft.c\" \(4203 characters\)
sed "s/^X//" >fft.c <<'END_OF_fft.c'
X/*
X *    fft.c
X *
X *    C Version 1.0 by Steve Sampson, Public Domain
X *
X *    This program is based on the work by W. D. Stanley
X *    and S. J. Peterson, Old Dominion University.
X *
X *    This program produces a Frequency Domain display
X *    from the Time Domain data input using the Fast Fourier Transform.
X *
X *    The REAL data is generated by the in-phase (I) channel and the
X *    IMAGINARY data is produced by the quadrature-phase (Q) channel of
X *    a Doppler Radar receiver.  The middle filter is zero Hz.  Closing
X *    targets are displayed to the right, and Opening targets to the left.
X *
X *    Note: With IMAGINARY data set to zero the output is a mirror image.
X *
X *    Usage:    fft  samples  input_data  output_data
X *    Where 'samples' is a power of two
X *
X *    Array Version for Turbo C 1.5
X */
X
X/* Includes */
X
X#include <stdlib.h>
X#include <stdio.h>
X#include <math.h>
X#include <alloc.h>
X
X/* Defines */
X
X#define    TWO_PI    ((double)2.0 * M_PI)
X
X/* Globals */
X
Xint    samples, power;
Xdouble    *real, *imag, max;
XFILE    *fpi, *fpo;
X
X/* Prototypes and forward declarations */
X
Xvoid    fft(void), max_amp(void), display(void);
Xint    permute(int);
Xdouble    magnitude(int);
X
X/* The program */
X
Xmain(argc, argv)
Xint    argc;
Xchar    *argv[];
X{
X    int    n;
X
X    if (argc != 4)  {
Xerr1:        fprintf(stderr, "Usage: fft samples input output\n");
X        fprintf(stderr, "Where samples is a power of 2\n");
X        exit(1);
X    }
X
X    samples = abs(atoi(argv[1]));
X    power =  log10((double)samples) / log10((double)2.0);
X
X    if ((real = (double *)malloc(samples * sizeof(double))) == NULL)
Xerr2:        fprintf(stderr, "Out of memory\n");
X
X    if ((imag = (double *)malloc(samples * sizeof(double))) == NULL)
X        goto err2;
X
X    if ((fpo = fopen(argv[3], "wb")) == (FILE *)NULL)
X        goto err1;
X
X    if ((fpi = fopen(argv[2], "rb")) != (FILE *)NULL)  {
X        fread(real, sizeof(double), samples, fpi);
X        fread(imag, sizeof(double), samples, fpi);
X        fclose(fpi);
X    }
X    else
X        goto err1;
X
X    fft();
X    max_amp();
X    display();
X
X    fwrite(real, sizeof(double), samples, fpo);
X    fwrite(imag, sizeof(double), samples, fpo);
X
X    fclose(fpo);
X}
X
X
Xvoid fft()
X{
X    unsigned i1, i2, i3, i4, y;
X    int     loop, loop1, loop2;
X    double     a1, a2, b1, b2, z1, z2, v;
X
X    /* Scale the data */
X
X    for (loop = 0; loop < samples; loop++)  {
X        real[loop] /= (double)samples;
X        imag[loop] /= (double)samples;
X    }
X
X    i1 = samples >> 1;
X    i2 = 1;
X    v = TWO_PI * ((double)1.0 / (double)samples);
X
X    for (loop = 0; loop < power; loop++)  {
X        i3 = 0;
X        i4 = i1;
X
X        for (loop1 = 0; loop1 < i2; loop1++)  {
X            y = permute(i3 / i1);
X            z1 =  cos(v * y);
X            z2 = -sin(v * y);
X
X            for (loop2 = i3; loop2 < i4; loop2++)  {
X                a1 = real[loop2];
X                a2 = imag[loop2];
X
X                b1 = z1*real[loop2+i1] - z2*imag[loop2+i1];
X                b2 = z2*real[loop2+i1] + z1*imag[loop2+i1];
X
X                real[loop2]      = a1 + b1;
X                imag[loop2]      = a2 + b2;
X
X                real[loop2 + i1] = a1 - b1;
X                imag[loop2 + i1] = a2 - b2;
X            }
X
X            i3 += (i1 << 1);
X            i4 += (i1 << 1);
X        }
X
X        i1 >>= 1;
X        i2 <<= 1;
X    }
X}
X
X/*
X *    Find maximum amplitude
X */
X
Xvoid max_amp()
X{
X    double    mag;
X    int    loop;
X
X    max = (double)0.0;
X    for (loop = 0; loop < samples; loop++)  {
X        if ((mag = magnitude(loop)) > max)
X            max = mag;
X    }
X}
X
X/*
X *    Display the frequency domain.
X *    The filters are aranged so that DC is in the middle filter.
X *    Thus -Doppler is on the left, +Doppler on the right.
X */
X
Xvoid display()
X{
X    int    c, n, x, loop;
X
X    n = samples / 2;
X
X    for (loop = n; loop < samples; loop++)  {
X        x = (int)(magnitude(loop) * (double)56.0 / max);
X        printf("%d\t|", loop - n);
X        c = 0;
X        while (++c <= x)
X            putchar('=');
X
X        putchar('\n');
X    }
X
X    for (loop = 0; loop < n; loop++)  {
X        x = (int)(magnitude(loop) * (double)56.0 / max);
X        printf("%d\t|", loop + n);
X        c = 0;
X        while (++c <= x)
X            putchar('=');
X
X        putchar('\n');
X    }
X}
X
X/*
X *    Calculate Power Magnitude
X */
X
Xdouble magnitude(n)
Xint    n;
X{
X    n = permute(n);
X    return (sqrt(real[n] * real[n] + imag[n] * imag[n]));
X}
X
X/*
X *    Bit reverse the number
X *
X *    Change 11100000b to 00000111b or vice-versa
X */
X
Xint permute(index)
Xint    index;
X{
X    int    n1, result, loop;
X
X    n1 = samples;
X    result = 0;
X
X    for (loop = 0; loop < power; loop++)  {
X        n1 >>= 1;            /* n1 / 2.0 */
X        if (index < n1)
X            continue;
X
X        result += (int) pow((double)2.0, (double)loop);
X        index -= n1;
X    }
X
X    return result;
X}
X
X/* EOF */
END_OF_fft.c
if test 4203 -ne `wc -c <fft.c`; then
    echo shar: \"fft.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0